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Sj^; ABSTRACT 

Q , We have developed a semi-empirical simulation of the extragalactic radio continuum sky suit- 

^ ■ able for aiding the design of next generation radio interferometers such as the Square Kilo- 

I meter Array (SKA). The emphasis is on modelling the large-scale cosmological distribution 

. of radio sources rather than the internal structure of individual galaxies. Here we provide a 

description of the simulation to accompany the online release of a catalogue of ~ 320 million 
simulated radio sources. The simulation covers a sky area of 20 x 20 deg^ - a plausible upper 
^ i limit to the instaneous field of view attainable with future (e.g. SKA) aperture array technolo- 

' gies - out to a cosmological redshift of z = 20, and down to flux density limits of 10 nJy 

at 151, 610 MHz, 1.4, 4.86 and 18 GHz. Five distinct source types are included: radio-quiet 
Tij" ' AGN, radio-loud AGN of the FRI and FRII structural classes, and star-forming galaxies, the 

CO ' latter split into populations of quiescent and starbursting galaxies. 

ly-^ , In our semi-empirical approach, the simulated sources are drawn from observed (or ex- 

■ trapolated) luminosity functions and grafted onto an underlying dark matter density field 

OO ! with biases which reflect their measured large-scale clustering. A numerical Press-Schechter- 

. style filtering of the density field is used to identify and populate clusters of galaxies. For 

economy of output, radio source structures are constructed from point source and elliptical 
sub-components, and for FRI and FRII sources an orientation-based unification and beaming 
model is used to partition flux between the core and extended lobes and hotspots. The ex- 
^ ' tensive simulation output gives users the flexibility to post-process the catalogues to achieve 

more complete agreement with observational data in the years ahead. The ultimate aim is for 
the 'idealised skies' generated by this simulation and associated post-processing to be fed to 
telescope simulators to optimise the design of the SKA itself. 

Key words: galaxies: active - galaxies:starburst - galaxies:luminosity function - galax- 
ies:radio continuum - cosmology:large-scale structure 
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1 INTRODUCTION 

The Square Kilometer Array (SKA) is the next generation radio 
telescope facility for the 21st century. Although it is still in the de- 
sign stage and unlikely to be fully operational before 2020, various 
pathfinder experiments will enter scientific service around 2010; 
one of these pathfinders will probably form the nucleus of the even- 
tual SKA and over the next decade grow in sensitivity from < 1 to 
100 per cent SKA. The scientific goals of the SKA have been de- 
scribed at length elsewhere (see e.g. Carilli & Rawlings 2004), but 
it is now clear that a new generation of more sophisticated science 
simulations is needed to optimise the design of the new telescopes 
and their observing programmes for the efficient realisation of these 



goals. With this in mind, we present here a new semi-empirical sim- 
ulation of the extragalactic radio continuum sky. The simulation is 
part of a suite of simulations developed under the European SKA 
Design Study (SKADS) initiative and referred to collectively as 
SKADS Simulated Skies (S^). It is one of two simulations aimed 
at simulating the extragalactic continuum and line-emitting radio 
sky, which together offer distinct yet complementary approaches to 
modelling the radio sky. The first approach, used by the simulations 
described in this paper, may be best described as 'semi-empirical', 
in the sense that the simulated sources are generated by sampling 
the observed (or extrapolated) radio continuum luminosity func- 
tions. The second approach, pursued by Obreschkow et al. (2008), 
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is dubbed 'semi-analytical' as it ascribes gas properties, star forma- 
tion and black hole accretion rates to galaxy haloes identified in an 
N-body simulation. Given their different emphases, the simulations 
are also commonly referred to as the 'continuum' (semi-empirical) 
and 'HI' (semi-analytical) simulations and, where appropriate, syn- 
ergies between them will be highlighted. Limitations in state-of- 
the-art N-body simulations mean that the semi-empirical simula- 
tions cover ^ 100 times more sky area than the semi-analytical 
ones. 

The purpose of the present paper is to describe the essential 
methodology and basic ingredients of the semi-empirical simula- 
tion (section 2); to demonstrate that its output satisfies some basic 
tests, such as reproducing the radio source counts, local luminos- 
ity function and angular clustering (section 3); to highlight some 
known deficiencies of the simulations (section 4). Simulation data 
products (source catalogues and images) can be accessed through a 
database web-interface (http://s-cube d.physics.ox.ac.uk ) where full 
details of data format can be found. 



2 DESCRIPTION OF THE SIMULATION 

There have been several previous attempts to model the radio 
source counts at the faint flux levels appropriate for the SKA and 
to produce simulated source catalogues, e.g. Hopkins et al. (2000), 
Windhorst (2003), Jackson (2004), Jarvis &. Rawlings (2004). Our 
new simulation builds on these efforts and, like them, is empirically 
based in the sense that it draws sources at random from the ob- 
served (or suitably extrapolated) radio continuum luminosity func- 
tions for the extragalactic populations of interest. In a major ad- 
vance over these works, however, we have attempted to model the 
spatial clustering of the sources instead of simply placing them on 
the sky with a uniform random distribution. This is an important 
feature because large scale structure measurements have emerged 
as a key SKA goal, stemming from the possibility to constrain the 
properties of dark energy through the detection of baryon acous- 
tic oscillations (BAG) (e.g. Abdalla & Rawlings 2005). The spatial 
clustering has been incorporated on large scales (3> 5/i~^ Mpc) by 
biasing the galaxy populations relative to an underlying dark mat- 
ter density field which is evolved under linear theory. Gn smaller 
scales, we have also used the underlying dark matter density field 
to identify and populate clusters of galaxies. 

Another major impetus for the development of the new sim- 
ulations has been the desire to generate an extensive suite of out- 
put catalogues which can be stored as a database and queried via 
a public web interface. The aim is to use such a tool to generate 
'idealised skies' which can then be fed to telescope simulators for 
the production and analysis of artifical datasets. Such a process is 
essential for optimising the design of the telescopes and their ob- 
serving programmes. 

2.1 Design requirements of the simulation 

The semi-empirical simulation has been designed with the key 
SKA science goals in galaxy evolution and large-scale structure in 
mind. The basic parameters are as follows: 

• Sky Area: 20 x 20 deg^, approximately matching the largest plau- 
sible instantaneous field of view of a future SKA aperture array. 
The latter will be used to perform HI redshift surveys out to 2 ~ 1.5 
over a significant fraction of the sky, in order to measure the galaxy 
power spectrum and the evolution of dark energy via baryon acous- 
tic oscillations (BAG) (see e.g. Abdalla & Rawlings 2005). 



• Maximum redshift: z = 20. Probing the "dark ages" in HI and 
radio continuum is another key SKA goal (see e.g. Carilli 2008), 
and the redshift limit of the simulation has accordingly been set to 
include the reionization epoch. 

• Flux density limits: 10 njy at the catalogue frequencies of 
151 MHz, 610 MHz, 1.4 GHz, 4.86 GHz and IB GHz, spanning 
the expected operational frequency range of the SKA. At 1.4 GHz, 
10 nJy is of the order of the expected sensitivity for a lOOh observa- 
tion with the nominal SKA sensitivity (A/Tsys = 20000 m^ K"\ 
350 MHz bandwidth). The flux limits are applied such that if any 
structural sub-component of a galaxy (i.e. core, lobes or hotspots 
in the case of a radio-loud AGN) has a flux density exceeding this 
limit at one or more of these frequencies, all sub-components of the 
galaxy are included in the output catalogue, regardless of whether 
they are all above the flux limit. In this paper, we neglect any in- 
trinsic or induced polarization properties of the radio emission, but 
we expect the first major revision of the catalogue to include such 
information. 

• Source populations: In addition to the classical double-lobed 
radio-loud AGN of the FRI and FRII classes, the simulation will 
also incorporate radio emission from 'radio-quiet' AGN and star- 
forming galaxies. The latter population comprises relatively quies- 
cent or 'normal' late-type galaxies, as well as more luminous and 
compact starburst galaxies. Together these two populations are ex- 
pected to dominate the radio source counts below a few 100 ^Jy at 
1.4 GHz. 

2.2 The theoretical framework 

The starting point for the simulation is a present-day {z = 0) 
dark matter density field, 5p/ p, defined on a cuboid grid of 
5/i"^Mpc comoving cells with the overall array having dimen- 
sions 550 X 550 X 1550 cells (the long axis defines the direc- 
tion of increasing redshift). An imaginary observer is situated at 
the centre of one face of the grid and looks out into the simulation 
volume over a solid angle AS7. The simulation consists in looping 
over all cells contained within AS7 (see illustration in Fig. [T]!. At 
the position of the ith cell, the comoving distance to the observer 
is inverted to yield an observed redshift, 21, for the assumed cos- 
mological model. For each source type (denoted by suffix j), the 
observed-frame spectral energy distribution is compared with the 
flux density limits of the simulation to derive a minimum lumi- 
nosity, -Lmin,j> for sources to be included in the output catalogues. 
The luminosity function is integrated down to this limit to yield 
the space density, (f){L > L]^^-^ ^, Zi). In the absence of any large- 
scale structure, the mean number of sources expected in the cell 
would be simply nj,xp j = (j){L > LJ„i,j j, Zi)AV^, where is 
the co-moving cell volume. To introduce large-scale structure into 
the distribution of the sources, we modify this as follows: 

nU,j = Ae''(^')^(^-)'''/''0(L > L^nj , ^i) AF, (I) 

where G(z) is the linear growth factor for fluctuations in the un- 
derlying dark matter density field, b(z) is the redshift-dependent 
bias of the population of interest, and yl is a normalisation factor 
which ensures agreement with the luminosity function when aver- 
aged over the largest scales. For a gaussian density field, it is easily 
shown that A = e"'''^^''^^^' " where cr^ is the cell variance 
of the density field. Simulated sources are generated by Poisson- 
sampling the luminosity function at L > ilninj with mean nj,xpj. 

The formalism of eqn. (I) is a somewhat adhoc, non-linear 
bias prescription which has the effect of amplifying the source den- 
sity in over-dense regions and depressing it in under-dense regions. 
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Figure 1. A 2-dimensional illustration of the basic geometry of the 
simulation. As discussed in section 2.2, the simulation consists in 
looping over all cells (each of which is5h^^ Mpc in size) contained 
within a cone of solid angle AQ, at the apex of which sits an imag- 
inary observer. At the position of the ith such cell (shaded), the 
co-moving distance to the observer is inverted to yield the redshift, 
2i, of the cell. 

When the exponent is small it reduces to a linear bias prescription, 
i.e. Sn/n ~ b{z)G{z)5p/ p. The use of a log-normal density field 
has the added benefit of accounting for a small amount of non- 
linear evolution in Sp/p and also prevents the unphysical situation 
p < (see Coles & Jones 1991). The choice of the 5h^^ Mpc cell 
size reflects a compromise between ensuring that it is: (i) large 
enough to keep 5p/p in the linear or quasi-linear regime; (ii) large 
enough to keep the total number of cells in the simulation volume 
within computer memory constraints; (iii) small enough to have 
sufficient mass resolution to identify a cluster of galaxies of mass 
W^*h~^ Mq (see section 2.8). The computation of the bias b{z) is 
addressed in section 2.7. 

2.3 The input cosmology 

We use a spatially flat cosmology with parameters Hq = 

70kmS"^ MpC"\ = 0.3, = 0.7, /baryon = 0.16, 

CT8 = 0.74 and 

^scalar 

= 1. The CAMB software (Code for 
Anisotropies in the Microwave Background; Lewis et al. 2000) is 
used to compute a linear isotropic power spectrum complete with 
BAO and an appropriate transfer function. The resulting z — 
power spectrum is sampled to generate Fourier modes of the den- 
sity field within the simulation volume, which is then fast Fourier- 
transformed to yield the z — dark matter density field, Sp/p, 
defined on the grid of 5h^^ Mpc cells. The latter is the starting 
point for the main simulation (section 2.2) and the cluster finding 
algorithm (section 2.8). 

2.4 Luminosity functions 

Here we describe the luminosity functions used in the simulation. 
For all populations, it has been unavoidably necessary to extrapo- 
late these in luminosity and redshift beyond the regimes in which 
they were observationally determined. We therefore expect that for 
many applications users will need to post-process the source cat- 
alogues to implement more realistic forms of high-redshift evo- 



lution, mostly by applying some form of negative space-density 
and/or luminosity evolution. Alongside the description of each lu- 
minosity function, we indicate what we consider to be the most rea- 
sonable form of high-redshift post-processing for each luminosity 
function, based on the observational data available at the present 
time. However, in order not to restrict future applications of the 
simulations as constraints at high redshift improve in the years 
ahead, the online database tool will give the user the choice be- 
tween applying either this default post-processing option or a user- 
defined form with space-density evolution and/or luminosity evo- 
lution decreasing in power-law or negative exponential form (or no 
post-processing at all). 

In all cases, the radio emission has been modelled without 
including polarization properties and without making any explicit 
allowance for the effects of the increasing energy density of the 
cosmic microwave background radiation (CMB) at high redshift. 
We have also made no attempt to model the Sunyaev-Zeldovich 
effect on either cluster or galactic-scales. 

• Radio-quiet AGN. Classical radio-loud AGN constitute only a 
small fraction of the overall AGN population; a more complete 
global census of AGN activity is provided by the hard X-ray AGN 
luminosity function (HXLF). Using the latter in conjunction with a 
radio:X-ray luminosity relation, Jarvis & Rawlings (2004) showed 
that radio-quiet AGN can make a significant contribution to the 
well-characterised upturn in the 1.4 GHz source counts below 
1 mjy, a feature which had previously been ascribed largely (al- 
though not without controversy) to star-forming galaxies (see e.g. 
Seymour et al. 2004). Recent multi-wavelength identification work 
on several deep fields appears to substantiate this claim (e.g. Simp- 
son et al. 2006; Seymour et al. 2008; Smolcic et al. 2008), suggest- 
ing that possibly as many as half the sources at the level of a few 
tens of p}y may be radio-quiet AGN. 

We therefore follow the Jarvis & Rawlings (2004) pre- 
scription for including radio-quiet AGN in our simulation. We 
use the Ueda et al. (2003) AGN HXLF for the instrinsic 
(i.e. absorption-corrected) 2-10 keV band and their model for 
luminosity-dependent density evolution (LDDE). Since the Ueda 
et al. HXLF does not include Compton-thick AGN (i.e. those with 
obscuring column densities A^h > 1-5 x 10^^ cm^'^ ), we increase 
the space density given by this luminosity function by 50 per cent 
(independent of X-ray luminosity) in a notional attempt to account 
for such sources. Although there is now compelling evidence from 
observations (e.g. Risaliti et al. 1999 at low redshift; Martinez- 
Sansigre et al. 2007 at high redshift) and X-ray background models 
(e.g. Gilli et al. 2007) that obscured Compton-thick AGN may be 
as abundant as unobscured AGN, the precise breakdown as a func- 
tion of intrinsic luminosity and redshift is not well-constrained. We 
therefore do not believe it worthwhile to attempt a more sophisti- 
cated correction for Compton-thick AGN at this stage and instead 
leave the matter open for refinement during post-processing. 

The intrinsic X-ray (2-10 keV; erg ) and radio luminosi- 
ties (1.4 GHz; W/Hz/sr) are assumed to be related as follows, as 
inferred from the correlation found by Brinkmann et al. (2000) for 
radio-quiet quasars: 

log(L2-10kev) = 1.012 log(Li.4GHz) + 21.3 (2) 

The correlation of Brinkmann et al. exhibits scatter of ~ ±1 dex 
at a given X-ray luminosity (and is further biased by the fact that 
it refers to a sample detected in both soft X-rays and at 1.4 GHz), 
but the extent to which this scatter is intrinsic to the physics of the 
emission process is not clear. It may be largely due to observational 
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errors, coupled with uncertain corrections for X-ray absorption and 
the effects of Doppler de-boosting of radio-core emission. How- 
ever, the fact that the slope of eqn. (2) is so close to unity strongly 
argues that this radio emission is a direct tracer of the nuclear activ- 
ity, and we therefore follow Jarvis & Rawlings (2004) in assuming 
a 1 : 1 relation with no scatter. 

The adopted luminosity limits of log(Li.4 ghz)=18.7 and 
25.7 [W/Hz/sr] span the range from low-luminosity AGN to the 
most luminous quasars. The sources are modelled as point sources 
with power-law radio spectra of the form oc v^"''' (see e.g. 
Kukula et al. 1998). The Ueda et al. dataset constrains the HXLF 
only out to z ~ 3 and assumes that the space density of all sources 
declines as (1 -|- z)~^'^ beyond z — 1.9, and this is assumed 
for the simulation. However, recent constraints from the Chandra 
Mult-wavelength project (CHAMP) (Silverman et al. 2007) have 
pushed the high-redshift constraints out to 2: ~ 5 and suggest 
that the decline in space density above z = 1.9 is steeper: 
(1 + Our default post-processing option is therefore to 

reduce the space density by a factor ( 1+z) "^'^^ above this redshift. 

• Radio-loud AGN. For the classical double-lobed radio-loud AGN, 
we use the 151 MHz luminosity function derived by Willott et 
al. (2001) from a compilation of low-frequency-selected samples. 
In particular, we use their 'model C luminosity function (adapted 
to our chosen cosmology) which consists of low and high luminos- 
ity components with different functional forms and redshift evo- 
lutions. The low-luminosity component consists of a Schechter 
function with an exponential cut-off above log L151 mhz — 
26.1 [W/Hz/sr] (in the = 1 cosmology of Willott et al.) and 
space density evolution of the form (1 + z)*"^ out to z = 0.706 
and constant thereafter. Willott et al. were not able to place con- 
straints on any high-redshift decline in this low-luminosity compo- 
nent of the luminosity function and our simulation therefore does 
not incorporate one. As a default post-processing option, however, 
we suggest a decline in the space density as (1 + z)^"^'^ above 
z = 2.5, matching that inferred by Jarvis & Rawlings (2000) for 
flat-spectrum quasars. 

The high luminosity component consists of Schechter func- 
tion with an inverted exponential cut-off at low luminosity, giving 
rise to a function which peaks at log L151 mhz — 27 [W/Hz/sr] (in 
the Willott et al. cosmology). The redshift evolution is modelled as 
a Gaussian peaking at 2; = 1.91 and no further post-processing is 
proposed or deemed necessary. These properties of the luminosity 
function make it natural to identify the low and high luminosity 
components crudely with FRI and FRII radio sources, respectively. 
We note, however, that Willott et al. pointed out that the luminosity 
at which the low- and high-luminosity components contribute 
equally to the luminosity function tends to be about 1 dex higher 
than the traditional FRI/FRII break luminosity. The adopted limits 
of integration for the luminosity function are: log L151 mhz = 20 
and 28 [W/Hz/sr] for the FRIs, and log L151 mhz = 25.5 and 
30.5 [W/Hz/sr] for the FRIIs (all in the simulation cosmology). 

• Star-forming galaxies. For star-forming galaxies, we essentially 
use the 1.4 GHz luminosity function derived from the IRAS 2Jy 
sample by Yun, Reddy & Condon (2001). This is modelled as a 
sum of two Schechter functions, with the high luminosity compo- 
nent dominant at log Li,4 ghz > 22 [W/Hz/sr], which equates to 
a 60pim luminosity of approximately 10^^ L©. For the purposes 
of assigning morphologies and SEDs, we follow Yun, Reddy & 
Condon and identify these low and high luminosity components 
with normal/quiescent late-type galaxies and starburst galaxies, re- 



spectively. This separation is borne out by the analysis of Takeuchi 
et al. (2003), who remeasured the IRAS 60/im luminosity func- 
tion and showed that these low and high luminosity components 
correspond to the cool and warm subpopulations (defined to have 
(3 = S'(100^m)/S'(60Atm) flux ratios (3 > 2.1 and < 2.1) respec- 
tively. 

We employ lower and upper luminosity limits of log 
L1.4 GHz = 17 and 25.5 [W/Hz/sr] for the normal galaxies, and log 
L1.4 GHz ~ 21 and 27 [W/Hz/sr] for the starbursts. The behaviour 
of the faint end of the radio luminosity function for the star-forming 
galaxies is a controversial but critical issue. The observational 
data of Yun, Reddy & Condon (2001) clearly diverge below their 
Schechter function fit below log L1.4 ghz = 19.3 [W/Hz/sr], an ef- 
fect the authors ascribe to measurement incompleteness. The star- 
forming galaxy luminosity function of Sadler et al. (2002), based 
on cross-matching the 1.4 GHz NVSS with the 2dFGRS, also peaks 
around log L1.4 ghz = 19 [W/Hz/sr] and may decline towards 
lower luminosities. However, other determinations of the luminos- 
ity function with higher degrees of faint end completeness (e.g. 
Mauch & Sadler 2007, Condon et al. 2002) show instead that the lu- 
minosity function just flattens below log Li,4 ghz ^ 20 [W/Hz/sr], 
with no evidence for an actual downturn. To take account of this, 
we simply assume that the z — Yun, Reddy & Condon lu- 
minosity function is flat below log L1.4 ghz ~ 19.6 [W/Hz/sr]. 
It has been shown that the local 1.4 GHz luminosity function 
of normal galaxies is (when suitably transformed) indistinguish- 
able from the luminosity function of 60/im-selected galaxies (e.g. 
Condon 1989,1992). We therefore assume that the radio emission 
serves as a reliable tracer of all star formation (albeit subject to the 
uncertainty in the faint end of the radio luminosity function), and do 
not follow the example of Hopkins et al. (2000)/Windhorst (2003), 
who supplemented the warm IRAS luminosity function in their sky 
simulations with an additional normal galaxy population derived 
from an optical luminosity function. 

For both components of the luminosity function we assume a 
form of cosmological evolution which, in a Universe with Qm ~ 1 
and Qa ~ takes the form of pure luminosity evolution of the form 
(1 + z)^'^ out to 2: = 1.5, with no further evolution thereafter. This 
is adapted to the cosmology of interest by shifting the luminosities 
and scaling the space density inversely with the ratio of the differ- 
ential co-moving volumes, dV/dz, for the two cosmologies. Evolu- 
tion of essentially this form (deduced for an Qm ~ 1 Universe) was 
inferred from studies of IRAS-selected galaxies (Rowan-Robinson 
et al. 1993) and for sub-mm-selected galaxies (Blain et al. 1999). 
More recent determinations of the evolution of the radio luminosity 
function of star-forming galaxies are consistent with this, mostly 
finding luminosity evolution of the form (1 + 2)*^ (Q = 2 — 3) 
and negligible density evolution (e.g. Seymour et al. 2004; Hop- 
kins 2004). Spitzer results on the evolution of the infrared (8- 
lOOO/xm) and far-infrared luminosity functions reveal a degener- 
acy between luminosity and density evolution which is consis- 
tent with our adopted form of evolution (Le Floc'h et al. 2005; 
Huynh et al. 2005), although constraints on the population with 
LiR < 10" L0 (equivalent to log I/i.4GHz = 22 [W/Hz/sr]) are 
limited at 2 > 0.5. 

More complex models of the evolution of star-forming galax- 
ies have been put forward to interpret mid- and far-infrared obser- 
vations (e.g. Franceshini et al. 2001; Lagache et al. 2004). Many 
such models break the population down into contributions from 
weakly or non-evolving normal galaxies, together with a popula- 
tion of starburst galaxies subject to strong evolution in space den- 
sity and luminosity. The model of Franceshini et al. (2001) is one 
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such model, which takes as its starting point a 12/im local luminos- 
ity function. However, when translated to 60^im with the assumed 
infrared SEDs, the resulting breakdown into normal and starburst 
galaxies (Fig. 1 1 of Franceschini et al. 2001) is completely different 
from the one we assume, with the starbursts comprising a decreas- 
ing fraction of the whole for Leo > 10^" L0. When translated to 
24^m (Gruppioni et al. 2005), the Franceschini et al. model is at 
variance with the redshift breakdown of Spitzer source counts re- 
ported by Le Floc'h et al. (2005). The Lagache et al. (2004) model 
decomposition of the local 60^m luminosity function into normal 
and starburst galaxies is similar to ours, but again there are discrep- 
ancies with respect to the Spitzer luminosity functions reported by 
Le Floc'h et al. (2005). In the absence of a definitive alternative 
model and to avoid the uncertainties arising from K-corrections in 
the mid-infrared, we prefer to base our predictions starting from the 
local 1.4 GHz luminosity function plus pure luminosity evolution 
for both sub-populations, as already discussed . 

The simulation was performed by assuming no further evolu- 
tion in the luminosity function beyond redshift z = 1.5. Taken at 
face value this is clearly unrealistic, but it gives the user freedom 
to impose their own form of high redshift decline by selectively 
filtering the catalogue during post-processing. At the present time, 
observations show that the star formation rate density is essentially 
constant from z ~ 1.5 to 2 ~ 4 and then to drops off sharply 
above this redshift (see e.g. Hopkins & Beacom 2006). As our 
default post-processing option, we model this high-redshift fall-off 
using the piece-wise power-law model of Hopkins & Beacom, in 
which the star-formation rate density falls off as (1 + above 
z = 4.8. 

• AGN/star-forming hybrid galaxies and double counting. The 
interplay between bursts of star formation and black hole accretion 
is of considerable interest for studies of galaxy evolution, and one 
aspect of this concerns hybrid galaxies in which both processes 
make significant contributions to the observed radio emission. 
Unfortunately, such objects cannot be isolated in the simulations 
due to the use of separate luminosity functions for the AGN and 
the starbursts, but they are almost certainly accounted for by virtue 
of the way in which the luminosity functions were constructed. A 
fraction of the luminous starbursts are likely to be mis-classified 
obscured AGN (or at least have a significant AGN contribution), 
and similarly some of the radio emission in the low luminosity 
AGN (especially the radio-quiet population) is likely to be related 
to star formation. The most likely upshot of this is that, in the 
simulation as a whole, the latter populations may have been double 
counted to a certain extent, thereby slightly overproducing the 
source counts. The semi-analytical simulations, in contrast, will 
include hybrid galaxies. A further instance of potential double 
counting concerns the fact that the hard X-ray luminosity function 
used for the radio-quiet AGN implicitly includes some contribution 
from the radio-loud AGN, so that the latter are in effect counted 
twice. However, this effect is likely dwarfed by the uncertainty 
in the Compton-thick correction factor applied to the X-ray 
luminosity function and we therefore make no further allowance 
for it. 



2.5 Radio-loud AGN: Unification, beaming, morphologies 
and radio spectra 

As described in section 2.4, radio-loud AGN are initially drawn 
from a 151 MHz luminosity function for steep-spectrum, lobe- 



dominated sources. We use this as the input parent population 
for an orientation-based unification model which enables us to 
assign source structures and radio spectra in a physically-motivated 
manner. Some of this procedure follows the prescription of Jackson 
& Wall (1999). The steps in our process are as follows: 

(i) Sources are assigned a true linear size, Dtruc, drawn at 
random from a uniform distribution [0,-Do(l + z)^^'"^}, where 
Do — 1 Mpc; this assumes that the sources expand with uniform 
velocity until they reach a size equal to the redshift-dependent 
upper envelope of the projected linear size distribution measured 
by Blundell, Rawlings & Willott (1999). 

(ii) The angle between the jet axis and the observer's line-of-sight 
is drawn from a uniform distribution in cos 6, and the jet axis is 
given a random position angle on the sky. 

(iii) The ratio of the intrinsic core to extended luminosities, 
defined at 1.4 GHz in the rest-frame, is given by Rch ~ 10^, 
where x is drawn from a Gaussian distribution of mean Xmcd and 
a = 0.5, truncated at abs(a;) > 10 to avoid numerical problems. 
The numerical values of Xmed we use are given at the end of this 
subsection. 

(iv) A relativistic beaming model is used to derive the ob- 
served coreiextended flux ratio, 7?obs = RchB{0), where 
B{9) = i[(l - /3cose)-2 + (1 + /3cos^)-2]; /3 = y(72-l)/7 
and 7 is the Lorentz factor of the jet. 

(v) The extended emission is modelled with a power-law spectrum 

oc u^"'"^'', while the core spectrum is modelled with some 
curvature: logFi, = aco + ad log v + ac2(log!^)^. For in 
GHz, Oci — 0.07, ac2 = —0.29, as measured by Jarvis & 
Rawlings (2000) from a sample of flat spectrum quasars (aco sets 
the normalization). The observational data on which these fits 
were based typically extend up to 10 or 20 GHz, implying that the 
model SEDs cannot be simply extrapolated to higher frequencies 
(e.g. to the WMAP bands above 20 GHz). 

(vi) For FRIs, the extended flux distribution on the sky is modelled 
as two coaxial elliptical lobes of uniform surface brightness, 
extending from the point source core, and each with a major axis 
length equal to half the projected linear size. The axial ratio of the 
lobes is drawn from a uniform distribution [0.2,1]. 

(vii) For FRIIs, the inner edges of the lobes are offset from the core 
by a distance / x PLS, where PLS is the projected linear size 
of the whole source and / is drawn from the uniform distribution 
[0.2,0.8]. A fraction /hs of the extended flux in the FRIIs is as- 
signed to point-source hotspots positioned at the ends of the lobes, 
where 

/hs = 0.4[log Li5i MHz - 25.5] ± 0.15, (3) 

in accordance with a correlation found by Jenkins & 
McEllin (1977) (the scatter is modelled with a uniform dis- 
tribution). The hotspots are assumed to have the same radio 
spectral shape as the rest of the lobe. 

(viii) In an attempt to model self-absorption within compact 
sources (the so-called GigaHertz Peaked Spectrum (GPS) sources) 
we apply a spectral turnover at a frequency, z^p, below which the 
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spectrum asymptotes to Fu oc v^'^ . We use the following relation, 
found by O'Dea (1998), between Vp and Dtruo: 

log !/p(GHz) = -0.21 - 0.65 log Dtruc(kpc). (4) 

The above formalism should be able to reproduce the observed 
variety of radio-loud AGN, including radio galaxies, steep and flat- 
spectrum spectrum radio quasars, and strongly beamed sources 
such as BL Lacs and blazars. The key parameters of the model 
are the core:lobe ratios (a^mcd) and jet Lorentz factors (7) for the 
FRI and FRII parent populations, and their possible dependence on 
ii5i MHz. There is a wealth of literature which could in princi- 
ple be used to constrain these parameters, but none of it is ideally 
suited to our needs due to a combination of selection effects and 
poorly defined asumptions. Instead, we have tuned the parameters 
to reproduce the observed multi-frequency radio source counts (as 
shown in section 3.2) at flux densities > 1 mjy where such sources 
dominate. We find [x^cd,l) = (—2.6,6) and (—2.8,8) for the 
FRI and FRIIs, respectively. For simplicity, we do not incorporate 
any luminosity dependence in these parameters. 

2.6 Star-forming galaxies: source sizes, radio spectra and HI 

As described in section 2.4, we split the star-forming galaxies into 
two sub-populations of 'normal' galaxies and 'starburst' galaxies, 
according to whether they are drawn from the low or high lumi- 
nosity Schechter function component of the Yun, Reddy & Con- 
don (2001) luminosity function, respectively. Here we discuss the 
assignment of radio spectra and structural properties to these pop- 
ulations. 

As reviewed by Condon (1992), the radio spectra of star- 
forming galaxies consist of two components: (i) thermal free-free 
emission with a luminosity directly proportional to the current pho- 
toionization rate and cc v^" with a ~ 0.1; (ii) non- thermal 
synchrotron emission from supemovae for which the power-law 
slope Qnt = 0.75 ± 0.1. The ratio of the total luminosity to that 
of the thermal component is not well constrained but is commonly 
assumed to be ~ 1 + 10(!/ghz)"'^ . In the presence of free- 
free absorption, acting on both the thermal and non-thermal compo- 
nents, the radiative-transfer solution for the overall radio spectrum 
takes the form: 

cx /.^[l - e-^«-][l + lOi^^'ir"]- (5) 

The free-free opacity, rg, scales with frequency as rg = 
{vs I vf''^ where vn is a normalisation parameter. For most normal 
galaxies = 0.003 — 0.01 GHz and for the observed frequencies 
of interest eqn. (5) can be recast as the sum of two power-laws: 

K'^GHz +10j^GSf- (6) 

For starburst galaxies, the more compact sizes and higher gas 
densities result in higher i^h values and we use — 1 GHz, as 
measured in the central regions of M82. For completeness, we add 
a far-infrared thermal dust component to the spectral energy distri- 
bution of star-forming galaxies using a modified black-body with 
T = 45K and P — 1.5. The dust SED is normalised using the 
far-infrared:radio correlation to relate the 60/^m and 1.4 GHz lumi- 
nosities (see Yun, Reddy & Condon 2001). This component has a 
discemable impact on the observed SED only at very high redshift 
z > 15 and at the highest output frequency (18 GHz). 

The assignment of disk sizes to the normal galaxies invokes a 
chain of reasoning which begins with the remarkably tight relation 
between HI mass (Mm) and disk size (Dm, in kpc, defined at an 



HI density of I M0 pc ^) measured by Broeils & Rhee (1997) for 
a sample of local spiral and irregular galaxies: 

log Mm = (1.96 ± 0.04) log Dhi + (6.52 ± 0.06). (7) 

We relate i\/Hi to Li.4ghz by combining the star formation 
rate-Li.4GHz relation of Sullivan et al. (2001) with the looser cor- 
relation between star formation rate and Mhi (Fig. 3, left panel in 
Doyle & Drinkwater 2006) to yield: 

log Mhi = 0.44 log Li.4GHz [W/Hz/sr] + 0.48 ± A, (8) 

where A is random scatter drawn from a normal distribution with 
a — 0.3. We derive a fiducial radio continuum disk diameter, 
L'cont, using the relationship given by Broeils & Rhee (1997) be- 
tween Dm and the optical absorption-corrected diameter, D^^ 
(measured at the 25th mag arcsec^'^ isophote): 

log Dcont = log Dm - 0.23 - log(l + z). (9) 

The (l + z) factor is not from Broeils & Rhee (1997) but an ap- 
proximate scaling of disk diameter with redshift suggested by cos- 
mological hydrodynamic simulations of disk formation (C. Power, 
private communication). For comparison, we note that Ferguson et 
al. (2004) suggest that galaxy radii scale as the inverse of the Hub- 
ble parameter, H{z)^^, as deduced from rest-frame UV observa- 
tions with the Hubble Space Telescope. 

For the starbursts, we assume the following continuum sizes: 

-Dcont = (l+2)^'kpC, (10) 

out to 2 = 1.5 and constant at 10 kpc for z > 1.5. This reflects the 
fact that local starbursts are compact (kpc-scale) whereas z ~ 2 
sub-mm galaxies are an order of magnitude larger in scale (see 
e.g. MERLIN interferometry, Chapman et al. 2004). Equation (8) is 
again used to assign a nominal HI mass to the starbursts, although 
the relation is unlikely to hold for such systems. In the output cat- 
alogues, star-forming disks are placed on the sky with random ori- 
entations. It should be noted that galaxies are catalogued in the out- 
put only according to their spatially integrated fluxes, and not on 
the basis of surface brightness. Users may thus model the source 
sizes and intensity profiles in a more sophisticated way during post- 
processing, should they so wish. 

2.7 Large-scale clustering and biasing 

The bias b[z) which appears in eqn. (I) is computed separately for 
each galaxy population using the formalism of Mo & White (1996). 
We make the simple assumption of assigning each population 
an effective dark matter halo mass which reflects its large-scale 
clustering. Our simulation does not have the mass resolution to 
directly resolve galaxy and group-sized haloes, so we are not 
literally assigning galaxies to individual haloes; our approach is 
thus distinct from the technique based on the 'halo occupation 
distribution function' which is used to populate dark matter haloes 
in N-body simulations, e.g. Benson et al. (2000). Whilst the as- 
sumption of a fixed halo mass for an entire population is reasonable 
at low redshifts, the formalism breaks down towards the higher 
redshifts of the simulation as b[z) increases, leading to a potential 
blow-up in the exponent of eqn. (I). To circumvent the excessively 
strong clustering which would otherwise result, the bias for 
each population is held constant beyond a certain cut-off red- 
shift (zcut). Halo masses are assigned to the populations as follows: 

• Radio-quiet AGN. The clustering of quasars and its redshift 
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evolution have been measured by the 2dFQZ survey (e.g. Croom 
et al. 2005) and found to be well described by a constant halo mass 
of 3 X lO^^/i^^ M© out to redshift hmits of the 2dFQZ {z = 2.5). 
We thus adopt this mass and set Zcut = 3. 

• Radio-loud AGN. The most extensive information on the 
cosmological clustering of radio sources has been derived from 
measurements of the angular clustering in the NVSS and FIRST 
radio surveys (see e.g. Overzier et al. 2003 and references therein). 
With assumptions about the redshift distribution for the radio 
sources, these clustering measurements can be inverted to yield 
3-D correlation lengths and dark halo masses. We use the discus- 
sion in Overzier et al. (2003) to fix our halo masses. For FRIs, 
we take Afhaio ~ 10^'^h~^ M©, which reflects the clustering 
of low-luminosity radio sources and L ^ L* E-type galaxies. 
For FRIIs, we take Afhaio = 10^'*/!.^^ M©, which reproduces 
the clustering of the powerful radio sources and matches the 
clustering evolution of z=0 L > L* E-type galaxies, z ~ 1 EROs 
and z ~ 0.55 LRGs. These masses are comparable to a small 
group and a small cluster, respectively, and are consistent with 
the findings of other studies which show that radio sources tend 
to reside in overdense, highly-biased environments. Since the 
resulting b{z) increases very rapidly with redshift, which would 
lead to an excessively inhomogenous clustering if left unchecked, 
we impose Zcut = 1-5. This is also the effective cut-off in the 
redshift distribution of clusters of galaxies (see next subsection), 
and the redshift beyond which, in any case, existing radio surveys 
do not constrain the clustering. 

• Star-forming galaxies. For the normal galaxies, we take 
Afhaio = lO^^/i~^M0 and Zcut = 3 which, as discussed by 
Overzier et al. (2003), reproduces the z=0 clustering of IRAS 
galaxies, L ~ 0.5 L* late-type galaxies and LBGs at 2: ~ 3. 
This population has a bias close to unity. For the starbursts, we 
take A/halo = 5 x 10^^ Mq, consistent with the claim of 
Swinbank et al. (2006) that submillimetre galaxies are progenitors 
of L > 3I/K E-type galaxies at z=0, and with the clustering 
measurements of Farrah et al. (2006) for mid-IR selected samples 
al z = 1.5 — 3. For such a group/cluster-sized mass, we again 
impose Zcut ~ 1.5. We note, however, that tuning the starburst 
clustering to match the high-redshift submillimetre population 
will not accurately reproduce the much weaker clustering of the 
local starburst galaxies, which is much lower than that of passive 
galaxies (e.g. Madgwick et al. 2003). 



2.8 Clusters of galaxies 

The simulation was designed for modelling large-scale structure 
on ^ 5h~^ Mpc scales. As discussed in section 2.7, the radio-loud 
AGN are heavily biased with respect to the underlying dark matter 
density field and are expected to lie preferentially in group/cluster 
environments. It was thus considered advantageous to be able to 
identify the clusters of galaxies in the simulation directly, since 
they are an important population in their own right for investiga- 
tions of cosmic magnetism and the Sunyeav-Zeldovich effect, for 
example. We attempted to do this whilst ensuring consistency with 
the framework for modelling large-scale structure outlined in sec- 
tion 2.2. 

This goal was accomplished through the use of a numerical 
Press-Schechter (Press & Schechter 1974) method to identify the 
cluster-sized haloes. The basic assumption is that any given cell is 



considered to be part of the largest mass halo which could have 
collapsed by the epoch under consideration. Under linear theory, a 
halo is defined to be collapsed when its extrapolated linear over- 
density reaches the value 5p/p = i5crit ~ 1.66.0 We thus filter the 
dark matter density field, 5p/p, on a range of successively smaller 
scales and at each filter step flag all pixels for which the smoothed 
density field exceeds (5crit /G{z) (where G{z) is as before the linear 
growth factor for the cell of interest, since we are once again look- 
ing out into the light cone). The next step is to locate the islands 
of inter-connected overdense cells and identify them as collapsed 
haloes, with a mass equal to that contained within the filter volume 
for the mean density of the Universe. The process is then repeated 
at the next lower mass scale, with the proviso that any cell already 
part of a larger mass structure is ignored. 

The finite cell size of the simulation introduces a discrete- 
ness in the filter size steps and hence in the masses of the result- 
ing clusters, although this becomes fractionally smaller at higher 
masses. The smallest mass filter includes seven 5/t~^ Mpc cells 
(corresponding to a cluster mass of 7.26 X 10^=^/1"^ Mq) and 
thereafter the filter grows in a quasi-spherical manner whilst re- 
taining symmetry along the xyz axes of the simulation volume. 
The resulting cluster masses are quantised to the following values: 
[0.726, 1.97, 2.80, 3.42, 5.91, 8.4] x 10"/i"^ Mq. This quantisa- 
tion of the cluster masses makes it very difficult to compare with 
continuous analytical cluster halo mass functions and observations 
of cluster counts; nevertheless, we attempted to do this by running 
the cluster-finder algorithm on a 10^/i^^Mpc^ volume at z = 0, 
assuming that the quantised masses are minimum estimates and that 
the true cluster masses actually lie in a range up to the next high- 
est mass level. The resulting mass function is shown in Fig. |2] for 
comparison with Press-Schechter (1974) and Sheth-Tormen (1999) 
mass functions; our mass function generally falls short of the lat- 
ter and some boosting of the masses will be necessary when post- 
processing the quantised cluster masses to generate a continuous 
mass distribution. As a crude example, we show that increasing all 
masses by a factor of 2 would bring the mass function into bet- 
ter agreement with the Sheth-Tormen function. A full investigation 
into the origin of this discrepancy is beyond the scope of this paper, 
and we simply note that it may result from a combination of sev- 
eral effects such as: (i) the use of a sharp-edged (top-hat) smoothing 
function for filtering the density field and the neglect of infall re- 
gions and inter-cell interpolation at the edge of the filter (e.g. the 
spherical filter enclosing a mass of 10^^ /i"^ Mq has a radius of 
14.2/i~^ Mpc, corresponding to just 2.8 cells, so systematic filter 
edge effects are likely to be significant); (ii) the fact that the ap- 
propriate mass density for the collapsed objects may slightly ex- 
ceed the assumed value (the mean density of the Universe); (iii) 
the reasons discussed by Sheth et al. (2001) for favouring an ellip- 
soidal collapse model - or, equivalently, a lower spherical collapse 
threshold - over the standard spherical Press-Schechter formalism. 

The areal density of clusters as a function of redshift is shown 
in Fig. [3] which can be compared in principle with observational 
results (e.g. van Breukelen et al. 2006) and N-body simulations 
(Evrard et al. 2002). However, the quantisation of the catalogued 

^ In the standard spherical collapse Press-Schechter model, (S^rit depends 
only weakly on the cosmological parameters - see e.g. Eke et al. (1996) - 
and is assumed to be independent of scale; Sheth et al. (2001) describe how 
an ellipsoidal collapse model with a scale- (i.e. mass) dependent threshold 
can improve the agreement with the mass function derived from N-body 
simulations, and that approximately the same results can be obtained in the 
spherical collapse model by reducing <5crit by — 16 per cent. 
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cluster masses again makes comparisons difficult. A proper com- 
parison will only be possible once conversion to a continous mass 
function has been carried out in post processing. Nevertheless, we 
note the results of van Breukelen et al. (2006), who report a sur- 
face density of ~ 10 clusters deg^^ for z — 0.5 — 1.5 and masses 
> 1O^''M0. By way of comparison, our total cluster catalogue 
shows a density of 5.6 clusters deg~^ above the 'boosted' mass of 
1.4 X 10"/i"^ Mq, 50 per cent of which are at z > 0.5. 

Having identified the clusters in the dark matter density field, 
we used the following procedure to populate them with galaxies. 
Any radio-loud or radio-quiet AGN in a cell within the smooth- 
ing filter radius of an identified cluster is assumed to be a mem- 
ber galaxy of that cluster. Star-forming galaxies of both types are 
placed without regard for the presence of clusters, i.e. just as de- 
scribed in section 2.2. We did not attempt to model the phase space 
distribution of cluster galaxies in a detailed manner; rather, galaxies 
were assigned a random orientation with respect to the cluster cen- 
tre and the physical radius was drawn at random from the uniform 
distribution [0, 0.5rviriai], where rviriai is the cluster virial radius. 
The latter was calculated from the cluster mass, Mci,i4 (in units of 
10^^ Mq), on the common assumption that the mean density 
within the cluster is ~ 100 times the background density of the 
universe (= ?,HIQ.m{1 + zf /^-kG) (see e.g. Eke et al. 1996): 

r^ai = lAMlll^{l + z)-^h-^Mpc (11) 

This is equivalent to assuming a 3-D galaxy density profile of 
the form p oc r^'^ within rviriai- Galaxies are offset in velocity from 
the cluster redshift by a random velocity drawn from a Gaussian 
distribution with standard deviation 0.8(7viriai- The calculation of 
the latter draws on the virial theorem ((Jviriai = GAfci/^viriai) to 
yield: 

CTviriai = 557M'/f4(l + zf^\m/s (12) 

Cluster membership is indicated in the output galaxy cata- 
logues with reference to a separate table listing the cluster prop- 
erties, giving users the freedom to populate the clusters in a more 
sophisticated manner during post-processing, should they so wish. 
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Figure 2. The crosses show the simulated cluster mass function 
derived from the 'quantised' cluster masses in a 10^/i~-^Mpc^ 
volume at z = 0. The solid and dotted lines show the Press- 
Schechter (1974) and Sheth-Tormen (1999) mass functions, respec- 
tively. It is clear that some additional boosting of the masses is nec- 
essary in post-processing when converting the quantised masses 
into a continuous distribution. The open triangles show the effect 
of artificially boosting all masses by a factor of 2, as discussed in 
section 2.8. 
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2.9 Redshift space distortions 

Within the linear-theory framework that we adopt, it is also possible 
to include the effect of peculiar velocities of galaxies on their mea- 
sured redshifts. Starting from the initial density field with Fourier 
modes (5k, the Fourier modes of the velocity field are given by: 



5vh{z) = - 



iH{z)5k{z) 
kG{z) 



AG{z) 



Az 



k 



(13) 



(Peebles 1976, Kaiser 1987, Peacock 1999). Each component of 
the peculiar velocity vector field is obtained by taking the Fourier 
transform of the vector modes in eqn. (13). The redshift distortion 
is then calculated at every point in the grid by projecting the pe- 
culiar velocity onto the line of sight and perturbing the original 
Hubble-flow redshift accordingly. When downloading source cat- 
alogues from the online database, users will have the choice be- 
tween selecting a pure Hubble-flow redshift, or a redshift which 
also includes the induced peculiar motions, the latter being the de- 
fault option. Note that this formalism does not include the effect 
on the peculiar velocities of non-linear growth at late cosmological 
epochs, but it does include the linear-theory distortion of measured 
clustering discussed by Kaiser (1987). 
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Figure 3. Surface densities of galaxy clusters lying at redshift z or 
higher, exceeding various mass thresholds. From top to bottom, the 
three lines correspond to mass thresholds of 7.0 x 10^^, 2.0 x 10^"*, 
and 3.4 x 10^^ h'^ Mq, respectively. When the 'quantised' masses 
are artificially boosted by a factor of 2 as described in Fig. |2l the 
same three lines now nominally correspond to mass limits 1.4 x 
10", 4.0 x 10" and 6.8 x 10"/i"^ Mq. 



3 SIMULATION OUTPUT AND BASIC TESTS 

Here we describe the output of the simulation and present a few 
consistency checks to demonstrate that it is in satisfactory agree- 
ment with existing observational constraints and can thus be ex- 
trapolated safely beyond this regime. The tests were performed in- 
dependently of the main simulation to provide maximum exposure 
to potential errors, using the raw simulation output (i.e. with no 
post-processing applied). 
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3.1 Simulation output : source catalogues 

The output of the simulation consists of two inter-linked 
catalogues, the radio source catalogue and the cluster cat- 
alogue, which can be accessed through a web-interface at 
[http://s-cubed.physics.ox.ac.uk Their structures are shown in Ta- 
bles 1 and 2, respectively. Each entry (row) of the radio source cat- 
alogue corresponds to a single structure (identified by the source in- 
dex) and lists the properties of that structure. Each galaxy is made 
up of one or more such structures, each having the same galaxy 
index, but different source indices, e.g. for each FRI radio source 
there are three structural components (point source core and two el- 
liptical lobes); FRIIs have five components (point source core, two 
elliptical lobes and two point-source hotspots). Radio-quiet quasars 
and star-forming galaxies all consist of a single structual compo- 
nent. A non-zero cluster index indicates that the galaxy is part of 
a cluster of galaxies, and the physical properties of the relevant 
cluster are listed in the cluster catalogue (Table 2). Note that the 
radio source table structure permits the specification of full polar- 
ization information (i.e. IQUV flux densities), but only the inten- 
sity I is generated by this simulation. Polarization information (i.e. 
QUV flux densities) will be supplied later by other members of the 
SKADS consortium. 

The total numbers of sources generated by the simulation are 
shown in Table 3, broken down into the numbers of galaxies and 
sub-structures (the two columns differ only for the FRIs and FRIIs). 
Note that the figures quoted in this table exceed by a few per cent 
those within the canonical 20 x 20 deg'^ simulation area, and that 
users should take care only to use sources within the latter area. 
These unavoidable edge effects arise from the discrete cells used 
in the simulation and diminish in importance towards the higher 
redshifts. The total number of galaxy clusters in the output cata- 
logue is 2248, of which 2201 are within the canonical 20 x 20 deg^ 
simulation area. 



3.2 Multi-frequency source counts 

Differential source counts have been generated at the five observa- 
tional frequencies (151 MHz, 610 MHz, 1.4 GHz, 4.86 GHz and 
18 GHz) and are shown in Fig. |4] In all cases, they are in good 
agreement with the existing observational data. We reiterate that 
the input luminosity functions for the simulations were inferred 
from observations at either 151 MHz or 1.4 GHz (depending on 
the galaxy type), so this agreement demonstrates that we are mod- 
elling the radio spectral energy distributions in a realistic fashion. 
For the radio-loud AGN, the beaming parameters were, however, 
derived in order to reproduce the source counts above 1 mjy 
(section 2.5). 

Concerning the faint end of the source counts which are cur- 
rently unconstrained by observational measurements, our model 
begins to diverge from the Hopkins et al. (2000)/Windhorst (2003) 
model below 1 ^Jy at 1.4 GHz and by 10 njy our counts are ap- 
proximately a factor of 10 below theirs. This is almost entirely due 
to the inclusion in their model of a 'normal galaxy' population de- 
rived from an optical luminosity function, supplementing the warm 
IRAS luminosity population. As Hopkins et al. point out, it is ques- 
tionable whether this additional population is needed as the normal 
galaxies may already be adequately represented by the tail of the 
IRAS luminosity function. As discussed in section 2.4, our model 
works on the latter assumption. There is, however, need for an addi- 
tional population of dwarf/irregular galaxies to account for the bulk 
of the HI mass function in range 10*^"** M© (section 3.7), but such 



galaxies are expected to be fainter in the radio continuum (com- 
pared with normal galaxies) and it is unclear how to incorporate 
them in these continuum simulations. 



3.3 Local luminosity functions 

The local 1.4 GHz luminosity functions for the individual galaxy 
types and for the radio source population as a whole have been con- 
structed and are shown in Fig. [5] 'Local' is defined here as meaning 
z < 0.3 and a lower flux density cut of 2.5 mJy at 1.4 GHz has 
been applied, in order to match the NVSS. The luminosity func- 
tions are in generally good agreement with those determined by 
Mauch & Sadler (2007) from a 6dF spectroscopic survey of NVSS 
radio sources. When the comparison is made for the AGN and star- 
forming galaxies separately, discrepancies appear at luminosities 
22 < logLi,4GHz < 23 [W/Hz/sr], which may be due to the ob- 
servational mis-classification of some of the AGN as star-forming 
galaxies, as discussed in section 2.4. 



3.4 Angular clustering 

As discussed in section 2, the incorporation of cosmological clus- 
tering in the galaxy distribution constitutes one of the main in- 
novations of our simulations. A simple test of this aspect is pro- 
vided by the angular correlation function, w{6), of radio sources, 
for which numerous measurements now exist from the FIRST and 
NVSS radio surveys at 1.4 GHz (see Overzier et al. 2003 and Blake 
& Wall 2002a,b). As discovered by the latter papers, the observed 
w{9) can be described by a double power-law: (i) a steep power law 
of the form w{d) — AG~^''^, dominant below 0.1 degrees, which 
is due to substructure (i.e. cores, lobes and hotspots) within indi- 
vidual radio galaxies; and (ii) a cosmological signal of the form 
w{9) = B6^"'^ which dominates on larger scales. According to 
Overzier et al. (2003), the amplitude of the cosmological clustering 
measured in NVSS/FIRST is essentially constant from 1.4 GHz 
flux density limits from 3-30 mJy, and then appears to increase by 
a factor of 10 above 200 mJy. However, a re-analysis of angular 
clustering in the NVSS, WENSS and SUMSS surveys by Blake 
et al. (2004) found no evidence that the amplitude depends on the 
1.4 GHz flux density over the range 10-200 mJy, although the as- 
sociated error bars are large. The 'w{9) measurements at these flux 
limits probe the clustering of just the powerful FRI/FRII sources 
and it would be necessary to extend them well below ImJy (or even 
100/xJy) in order to discern the impact of the star-forming galaxies 
on the clustering signal in an unidentified radio survey (Wilman et 
al. 2003). 

Our w{9) measurements are shown in Fig.|6]for 1.4 GHz flux 
limits of 3 and 10 mJy. Qualitatively, the simulations reproduce 
the observed double component structure in w{0) with the break 
in the observed position (0.1 degrees) and the correct slope and 
amplitude for the cosmological component. There are some devia- 
tions from a power-law form in the multi-component contribution 
to ■w{6) but these just reflect our simplistic modelling of the internal 
structures of individual radio sources. The amplitude of the simu- 
lated cosmological clustering signal at 3mJy is in agreement with 
the 'base level' amplitude measured by Overzier et al. (2003) at 
3-30 mJy. At higher flux densities, however, the simulated ampli- 
tude appears to increase rather faster than observed. This is likely 
to be the result of the relatively crude bias prescription in the simu- 
lations, in which a single halo mass is assigned to each population, 
independent of luminosity and redshift: in reality, there may be be 
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Figure 4. Source counts generated from the simulation catalogue at the five frequencies, broken down according to source type. The red lines 
are polynomial fits to the observed source counts from the following papers: 151 MHz, Willott et al. (2001); 610 MHz, Gam et al. (2008); 
1.4 GHz, Hopkins et al. (2003); 4.86 GHz, Gregory et al. (1994); in place of 18 GHz data we plot the 9C survey 15GHz source counts from 
Waldram et al. (2003). In all cases, the simulated source counts are within the scatter of the observations. 
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Table 1. Structure of the radio source catalogue 



Column 


Fortran Format 


Content description 




Tin 


oource inuex 


2 


16 


f~'liictpr itiHpY 


3 


110 


Galaxy index 


4 


14 


^tnr-TrirmatiriTi tvnf inHpv — nr\ 1 — nnrmal onlnw '7 — stnrhiirct"! 


5 


14 


AGN-fvnp index fO-no AGN l-ROO 7-FRT 3-FRTT 4-GPS1 


g 


14 


Ntniptiirp-tvnp ^^ — pnrp v — Irihp — nntcnrvt A, — Hict^ 


7 


F11.5 


RA {degrees) 


8 


F11.5 


DEC (degrees) 


9 


F10.3 


Comoving distance (Mpc) 


10 


F10.6 


Redshift 


11 


F9.3 


Position angle (rad) for elliptical substructures 


12 


F9.3 


Major axis (arcsec) 


13 


F9.3 


Minor axis (arcsec) 


14-17 


F9.4 


log(I), log(Q), log(U), log(V) flux densities (Jy) @ 151 MHz 


18-21 


F9.4 


log(I), log(Q), log(U), log(V) flux densities (Jy) @ 610 MHz 


22-25 


F9.4 


log(I), log(Q), log(U), log(V) flux densities (Jy) @ 1400 MHz 


26-29 


F9.4 


log(I), log(Q), log(U), log(V) flux densities (Jy) @ 4860 MHz 


30-33 


F9.4 


log(I), log(Q), log(U), log(V) flux densities (Jy) @ 18000 MHz 


34 


F9.4 


log (Mjji/Mq) (for star-forming galaxies only) 


35 


F9.4 


cos (viewing angle) (relative to jet axis; FRI and FRIl only) 



Table 2. Structure of the cluster catalogue 



Column 


Format 


Content description 


1 


Long integer 


Cluster index 


2 


Float 


RA (degrees) 


3 


Float 


DEC (degrees) 


4 


Float 


Redshift 


5 


Float 


Cluster mass {h^^ Mq) 


6 


Float 


Cluster virial radius (/i~^ Mpc) 


7 


Float 


Cluster velocity dispersion ( km s"'^ ) 



observational data from the 3CRR, 6CE and 7CRS surveys (e.g. 
Blundell, Rawlings & Willott 1999; BRW). Flux density cuts have 
been applied to the simulated sources to match these surveys. De- 
spite the simplicity of the algorithm for applying true linear sizes 
to our simulated sources (section 2.5), there is good agreement be- 
tween the projected linear sizes of the simulated powerful radio 
sources and those in real low-frequency-selected redshift surveys. 
In both plots of Fig. [8] the main differences between the simulated 
and real data points are readily explicable as being due either: (i) 
to the differences in area between the semi-empirical simulation 
and the real datasets; or (ii) to known inadequacies in the simula- 
tion technique. Examples of case-(i) differences include; the semi- 
empirical survey has a sky area a factor ~ 35 smaller than 3CRR 
(0.12 vs 4.239 sr), explaining why no simulated sources appear in 
the upper dex of radio power (i.e. 28 < log Lisi < 29 [W/Hz/sr]) 
in Fig. [8] compared with 37 3CRR sources; the small number of 
real objects beyond 2 ~ 3 in Fig. [8] (to which 3CRR is insensi- 
tive) reflects the fact that the simulation covers 5.5 times the area 
of the 7CRS (0.12 vs 0.022 sr). An example of a case-(ii) differ- 
ence is that the simulation hard-wires a Dtruc = 1 Mpc cut-off 
in size, which declines with redshift, whereas a few larger (giant) 
sources larger than this redshift-dependent limit are seen in the real 
data out to 2: ~ 2. There is clearly much scope for ameliorating 
such problems by applying more sophisticated radio sources evo- 
lution models, although as emphasised by BRW this needs to be 
done in a self-consistent way that simultaneously models the space- 



a smooth transition in halo mass from the FRIs (for which we as- 
sumed 10"/i-^ Mq) to the FRIIs (modelled with 10"/t~^ Mq). 

3.5 Spectral index distributions 

The distributions of observed radio spectral index between 
610 MHz and 1.4 GHz are shown in Fig. [7] broken down into 
the separate distributions for the star-forming galaxies and the 
AGN (spectral index a is here defined with the convention oc 
i/^"). Flux density limits of Sgiomhz > 360/iJy and S'i.4GHz > 
200^Jy were applied to match the selection criteria of Garn et 
al. (2008), whose spectral index distribution is shown for compar- 
ison. There are clearly discrepancies between the two, particularly 
the failure of the simulations to reproduce the flat spectrum tail of 
the observed distribution, which might be partially due to measure- 
ment errors at the flux density limit in the observations and also in 
part due to an excessive amount of high-frequency curvature in the 
assumed spectra of the radio-loud AGN cores (as defined in sec- 
tion 2.5). Fig. 7 also shows the effect of gaussian-smoothing the 
simulated spectral index distribution during post-processing. 

3.6 P-D and D-z diagrams 

In Fig. [8] we compare the 151 MHz radio luminosity-redshift and 
projected linear size-redshift relations of the simulated sources with 
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Table 3. Numbers of galaxies and their structural sub-component sources in the output catalogues of the 20 X 20 deg-^ semi-empirical simulations 



Source type Number of galaxies (10 ) Structural components (10 ) 



Radio-quiet AGN 


36.1 


36.1 


FRI 


23.8 


71.4 


FRII 


0.00235 


0.012 


Normal galaxies 


207.8 


207.8 


Starburst galaxies 


7.26 


7.26 



and-time-dependent spectral indices of the extended radio-emitting 
components and also models variations in the radio source environ- 
ments. 



3.7 HI mass function 

The main focus of these simulations is the radio continuum emis- 
sion. As discussed in section 2.6, however, we can use the loose 
correlation between HI mass and star formation rate to assign HI 
masses to the star-forming galaxies. Although this relation was 
measured for the 'normal galaxy' population, we also apply it to the 
starburst galaxies, even though it is unlikely to hold for the most ex- 
treme objects. We did not attempt to assign HI masses to the AGN 
populations because there are no simple observationally-based pre- 
scriptions for doing so. Moreover, the AGN are most likely to re- 
side in early-type galaxies, whose contribution to the local HI mass 
function is significantly less than that of the late-types (Zwaan et 
al. 2003). 

The z = — 0.1 HI mass function of the star-forming galax- 
ies in the simulation is shown in Fig [9] where it is compared 
with that from the HI Parkes All Sky Survey (HIPASS) (Zwaan 
et al. 2003). The excess contribution in the simulations over the 
HIPASS fit above 10^°-^ M© is due to the contribution of the star- 
burst galaxies which, as discussed above, are likely to have lower 
HI masses than those implied by the HI mass-star formation rate 
relation of the normal galaxies. The divergence between the sim- 
ulations and the HIPASS fit at lower masses is chiefly due to the 
absence in our simulation of a dwarf/irregular galaxy population; 
Zwaan et al. showed that this population (morphological type "Sm- 
Irr") accounts for an increasing proportion of the HI mass function 
as the mass is reduced below 10^" M©, and essentially for all of 
it in the range 10^^** M©. However, such a population is proba- 
bly not accounted for in our normal galaxy radio luminosity func- 
tion, even though the lower integration limit for the latter of log 
Li.4 = 17 [W/Hz/sr] equates to a very modest star formation rate 
of around 10^'^ M© yr~^ . Compared with normal galaxies, dwarf 
galaxies are gas-rich objects, with higher ratios of HI mass to star 
formation rate (see e.g. Roberts & Haynes 1992), possibly because 
conditions in the dwarf galaxies do not meet disc instability criteria 
necessary for the onset of star formation. 

A much fuller treatment of HI is provided by the SKADS 
semi-analytic simulations of Obreschkow et al. (2008). 



4 DEFICIENCIES OF THE SIMULATION 

Here we reiterate what we regard as the most important weaknesses 
in our modelling approach that users should be aware of. Some of 
these issues can be circumvented with appropriate post-processing, 
whereas others are hard-wired into the simulation and therefore of 



a more fundamental nature. 

• The use of extrapolated luminosity functions. The extrapolation 
of luminosity functions beyond the regimes of luminosity and 
redshift in which they were determined is unavoidable when 
attempting to simulate a next generation facility with a quantum 
leap in sensitivity over existing telescopes. The most important 
aspect of this concerns the faint end of the normal galaxy radio 
luminosity function, which is controversial even in the regime 
already observed. As discussed in section 2.4, we assume that the 
luminosity function flattens below log Li.4ghz = 19.6 [W/Hz/sr], 
as determined by Mauch & Sadler (2007), whilst others (e.g. 
Hopkins et al. 2000) have supplemented this with an additional 
population derived from an optical luminosity function. We also 
assumed that the luminosity function of the star-forming galaxies 
does not evolve further in redshift from z = l.b out to the redshift 
limit, 2 = 20. This is clearly unrealistic, but it gives the user 
full freedom to implement any particular form of decline in the 
star formation rate as a post processing task, as multi-wavelength 
constraints on high-redshift star formation accrue in the coming 
years. For example, stronger forms of high-redshift decline in 
the space density can be implemented by random sampling the 
existing catalogue. For all luminosity functions, we have indicated 
the default post-processing option for negative high-redshift 
evolution, and on the web database users will have the choice 
between implementing this or some other form. 

• The lack of star-fonning/AGN hybrid galaxies and double 
counting. The use of separate luminosity functions for the AGN 
and star-forming galaxies is a fundamental design limitation of 
this simulation, and it prevents us from explicitly modelling hybrid 
galaxies where both processes contribute to the radio emission. 
However, such galaxies are implicitly present in the simulation 
for the reason that at the faint end of the AGN radio luminosity 
function, star formation may make a non-negligible contribution 
to the radio emission. Similarly, some of the powerful starbursts 
may contain sizeable AGN contributions. These effects will lead to 
a certain amount of double counting, perhaps over-producing the 
faint source counts. 

• The lack of small-scale non-linear clustering. Due to the use 
of 5h~^ Mpc cells and a linear-theory power spectrum, this 
simulation does not model small-scale non-linear clustering in 
a satisfactory manner. This is not considered a serious problem 
because the simulation was primarily designed to model the 
clustering on large-scales, which is adequately described by our 
linear theory prescriptions. The 'semi-analytical' simulations of 
Obreschkow et al. (2008), based on N-body simulations, provide a 
much better description of small-scale non-linear clustering. 

• The treatment of galaxy clusters. Clusters were identified in 



A semi-empirical simulation of the extragalactic radio continuum sky 13 



r\l 
I 




CD 
I 



20 22 24 26 

log^JL, , (Watts Hz"^ sr^)] 




19 20 21 22 23 24 

log,o[L, ^ (Watts Hz"' Sr')] 

Figure 5. The uppermost plot shows the contributions of the various 
populations to the local (2 < 0.3) 1.4 GHz luminosity function in 
the simulation, for flux densities > 2.5 mJy. The data points show 
the luminosity function derived for the same flux density limit for 
NVSS sources by Mauch & Sadler (2007). The middle and lower 
plots show the comparison for the AGN and star-forming galaxies 
separately, with the grey regions indicating the Poisson errors on 
the simulations. 
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Figure 6. The angular 2-point correlation function, w{0), measured 
in our simulations for 1.4 GHz flux limits of 3mJy and lOmJy. The 
lines overplotted in both cases are the cosmological (dashed line) 
and multi-component (dot-dashed line) contributions to w{0) mea- 
sured by Blake et al. 2004 for a lOmJy flux limit. The error bars are 
Poisson error bars, which are not strictly correct for clustered data. 



the mass density field using a numerical Press-Schechter filtering 
method. The coarse sampling of the density field and the sharpness 
of the real-space mass filter result in a quantisation of the cluster 
masses. Users with particular interest in the clusters would thus 
need to post-process the masses to obtain a continuous distribution 
to match the Sheth-Tormen (1999) mass function (which provides 
a better fit to N-body simulations at the high mass end than the 
standard Press-Schecther formalism). It should also be noted that 
galaxies were assigned to the clusters in a very simplistic fashion, 
with no account of infall regions, for example. A more fundamen- 
tal design limitation of the simulation is that the underlying galaxy 
luminosity functions are assumed to be the same in the cluster and 
field environments. 



5 CONCLUSIONS 

The semi-empirical simulation described in this paper has been 
designed to be used in an interactive fashion for optimising 
the design of the SKA and its observing programmes in or- 
der to achieve its scientific goals. For this reason, we have 
made the galaxy catalogues available via a web-based database 
(http://s-cubed.physics.ox.ac.uk;l, which users can interactively 
query to generate 'idealised radio skies' according to their require- 
ments. Such skies can then be fed to (software) telescope Simula- 
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Figure 7. The upper plot shows the observed 610 MHz-1.4 GHz 
spectral index distribution for our simulation, broken down into 
separate populations of AGN and star-forming galaxies (the spike 
around 0.7-0.8 is due to the radio-quiet AGN and the extended 
emission of the radio-loud AGN, which are modelled without scat- 
ter). The lower plot shows an observed distribution derived by Garn 
et al. (2008) and the effect of gaussian smoothing the simulated dis- 
tribution, in an attempt to better match the observations. 
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Figure 8. A comparison of the projected linear size-redshift and ra- 
dio power- projected linear size relations of the simulated sources 
(open points) with those in the 3CRR, 6CE and 7CRS surveys 
(filled points). Flux density cuts have been applied to the simulated 
sources to mimic these surveys. See text for discussion. 



tors to facilitate the development of data processing and calibration 
pipeline routines. The extent to which the scientific content of the 
simulations can then be recovered at the end of such a simulation 
chain will enable observing programmes to be optimised for vari- 
ous models of the telescope hardware. In parallel with such work, 
we also anticipate that the simulations will be a valuable tool in 
the interpretation of existing radio surveys with abundant follow- 
up data. 
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Figure 9. The solid line shows the HI mass function derived from 
the star-forming galaxies in the simulation over the redshift range 
2 = — 0.1. The dashed line is the local HIPASS mass function 
from Zwaan et al. (2003). 



REFERENCES 

Abdalla F.B., Rawlings S., 2005, MNRAS, 360, 27 

Benson A.J., Cole S., Frenk C.S., Baugh CM., Lacey C.G., 2000, MNRAS, 
311, 793 

Blain A., Small I., Ivison R.J., Kneib J.-R, 1999, MNRAS, 302, 632 
Blake C, Wall J., 2002a, MNRAS, 329, 37 



Blake C, Wall J., 2002b, MNRAS, 337, 993 
Blake C, Mauch T., Sadler E.M., 2004, MNRAS, 347, 787 
Blundell K.M., Rawlings S.G., Willott C.J., 1999, AJ, 117, 677 
Brinkmarm W., Laurent-Muehleisen S.A., Voges W., Siebert J., Becker 
R.H., Brotherton M.S.. White R.L., Gregg M.D., 2000, A&A, 356, 445 
Broeils A.H., Rhee M.-H.. 1997. A&A, 324, 877 
Carilli C, Rawlings S. (eds), 2004, NewAR, 48 



A semi-empirical simulation of the extragalactic radio continuum sky 15 



Carilli C, 2008, to appear in Proceedings of Science, 'From planets to dark 
energy: tlie modem radio Universe' ed. Beswick (astro-pli/0808.1727) 

Chapman S.C., Smai! I., Windliorst R., Muxlow T., Ivison R.J., 2004, ApJ, 
611, 732 

Coles R, Jones B., 1991, MNRAS, 248, 1 
Condon J.J., 1989, ApJ, 338, 13 
Condon J.J., 1992, ARAA, 30, 575 

Condon J.J., Cotton W.D., Broderick J.J., 2002, AJ, 124, 675 
Croom S.M., et al., 2005, MNRAS, 356, 415 
Doyle M.T., Drinkwater M.J., 2006, MNRAS, 372, 977 
Eisenstein D.J., Seo H.-J., White M., 2007, ApJ, 664, 660 
Eke V.R., Cole S.M., Frenk C.S., 1996, MNRAS, 282, 263 
Evrard A., et al., 2002, ApJ, 573, 7 
FaiTah D., et al., 2006, ApJ, 641, L17 
Ferguson H.C., et al., 2004, ApJ, 600, 107 

Franceschini A., Aussel H., Cesarsky C.J., Elbaz D., Fadda D., 2001, A&A, 
378, 1 

Gam T., Green D.A., Riley J.M., Alexander R, 2008, MNRAS, 383, 75 
GiUi R., Comastri A., Hasinger G., 2007, A&A, 463, 79 
Gregory RC, Vavasour J.D., Scott W.K., Condon J.J., 1994, ApJS, 90, 173 
Gruppioni C, Pozzi F, Lari C, Oliver S., Rodighiero G., 2005, ApJ, 618, 
LI 

Hopkins A.M.. Windhorst R., Cram L., Ekers R., 2000, ExA, 10, 419, 

( astro-pl-i/9906469 I 
Hopkins A.M., Afonso J., Chan B., Cram L.E., Georgakakis A., Mobasher 

B., 2003, AJ, 125, 465 
Hopkins A.M., 2004, ApJ, 615, 209 
Hopkins A., Beacom J.F, 2006, ApJ, 651, 142 
Jackson C.A., 2004, NewAR, 48, 1187 
Jackson C.A., Wall J.V., 1999, MNRAS, 304, 160 
Jarvis M.J., Rawlings S.G., 2000, MNRAS, 319, 121 
Jarvis M.J., Rawlings S.G., 2004, NewAR, 48, 1 173 
Jenkins CM., McEllin M., 1977, MNRAS, 180, 219 
Kaiser N.. 1987, MNRAS, 227, 1 

Kukula M.J., Dunlop J.S., Hughes D.H., Rawlings S., 1998, MNRAS, 297, 
366 

Lagache G., et al., 2004, ApJS, 154, 112 
Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473 
Martinez-Sansigre A., et al., 2007, MNRAS, 379, L6 
Mauch T., Sadler E.M.. 2007. MNRAS, 375, 931 
Madgwick D.S., et al., 2003, MNRAS, 344, 847 
Mo H.J., White S.D.M.. 1996. MNRAS, 282, 347 
Obreschkow D., et al., 2008, in preparation 
O'Dea C, 1998, PASR 110, 493 

Overzier R.A., Rottgering H.J.A., Rengelink R.B., Wilman R.J., 2003, 
A&A, 405, 530 

Peacock J. A., 1999, 'Cosmological Physics', publ. Cambridge University 
Press 

Peebles PJ.E., 1976, ApJ, 205, 318 

Press W.H., Schechter S., 1974, ApJ, 187, 425 

Risaliti G., Maiolino R., Salvati M., 1999, ApJ, 522, 157 

Roberts M.S., Haynes M.R, 1992, ARA&A, 32, 115 

Rowan-Robinson M., Benn C.R., Lawrence A., McMahon R.G., Broadhurst 

T.J., 1993, MNRAS, 263, 123 
Sadler E.M., et al., 2002, MNRAS, 329, 227 
Seymour N., McHardy I.M., Gunn K.F, 2004, MNRAS, 352, 131 
Seymour N., et al., 2008, to appear in MNRAS (astro-ph/0802.4105) 
Sheth R.K., Tormen G., 1999, MNRAS, 308, 1 19 
Sheth R.K., Mo H.J., Tormen G., 2001, MNRAS, 323, 1 
Smolcic v., et al., 2008, in press at ApJS (astro-ph/0803.0997) 
Silverman J.D., et al., 2007. submitted to ApJ (astro-ph/0710.2461) 
Simpson C, et al., 2006, MNRAS, 372, 741 

Sullivan M., Mobasher B., Chan B., Cram L.. Ellis R., Treyer M., Hopkins 

A.,2001, ApJ, 558, 72 
Swinbank A.M., Chapman S.C., Small I., Lindner C, Borys C, Blain A.W., 

Ivison R.J., Lewis G.F, 2006, MNRAS, 371, 465 
Takeuchi T.T., Yoshikawa K.. Ishii T.T., 2003, ApJ, 587, L89 
Ueda Y, Akiyama M., Ohta K., Miyaji T., 2003, ApJ, 598, 886 



van Breukelen C, et al., 2006, MNRAS, 373, L26 

Waldram E.M., Pooley G.G., Grainge K.J.B., Jones M.E., Saunders R.D.E., 

Scott RE., Taylor A.C., 2003, MNRAS, 342, 915 
Willott C.J., Rawlings S.G., Blundell K.M., Lacy M., Bales S.A., 2001, 

MNRAS, 322, 536 
Wilman R.J., Rottgering H.J.A., Overzier R.A., Jarvis M.J., 2003, MNRAS, 

339, 695 

Windhorst R., 2003, NewAR, 47, 357 

Yun M.S,. Reddy N.A., Condon J.J., 2001, ApJ, 554, 803 

Zwaan M.A., et al., 2003, AJ, 125, 2842 



